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This article is composed of two parts; In the first part (Sec. 
1), the ultra-large-scale electronic structure theory is reviewed 
for (i) its fundamental numerical algorithm and (ii) its role in 
nano-material science. The second part (Sec. 2) is devoted to 
the mathematical foundation of the large-scale electronic struc- 
ture theory and their numerical aspects. 

1 Large-scale electronic structure theory and 
nano-material science 

1 . 1 Overview 

Nowadays electronic structure theory gives a microscopic foun- 
dation of material science and provides atomistic simulations 
in which electrons are treated as wavefunction within quantum 
mechanics. An example is given in the upper left panel of Fig.l. 
For years, we have developed fundamental theory and program 
code for large-scale electronic structure calculations, particu- 
larly, for nano materials. [1-6] The code was applied to several 
nano materials with 10^-10^ atoms, whereas standard electronic 
structure calculations are carried out typically with 10^ atoms. 
Two application studies, for silicon and gold, are shown in the 
right panels and the lower left panels of Fig.l, respectively. Now 
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the code is being reorganized as a simulation package, named as 
ELSES (Extra-Large-Scale Electronic Structure calculation), for 
a wider range of users and applications in science and industry. 

[6] 




Figure 1: Upper left panel: Example of calculated electronic wavefunction 
on a silicon surface (A 'yr-type' electron state on Si(lll)-2xl surface, given by 
a standard electronic structure calculation). Right panels: Application of our 
code to fracture dynamics of silicon crystal. [2] In results, the fracture path 
is bent into experimentally-observed planes (right upper panels) and and 
reconstructed Si(lll)-2xl surfaces appear with step formation (right lower 
panels). Lower left panels: Application of our code to formation process 
of helical multishell gold nanowire [5] that was reported experimentally. [7] 
A non-helical structure is transformed into a helical one (left and middle 
panels). The section view (right panel) shows a multishell structure, called 
'11-4 structure', in which the outer and inner shells consist of eleven and four 
atoms, respectively. 
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1.2 Methodology 



Our methodologies contain several mathematical theories, as 
Krylov-subspace theories for large sparse matrices. Here we fo- 
cus on a solver method of shifted linear equations, called 'shifted 
conjugate-orthogonal conjugate- gradient (COCG) method' [3] [4]. 
A quantum mechanical calculation within our studies is conven- 
tionally reduced to an eigen- value problem with a real-symmetric 
NxN matrix H (Hamiltonian matrix), which will cost an 0{N^) 
computational time. In our method, instead, the problem is re- 
duced to a set of shifted linear equations; 

{z^^^I - H) x{z^^^) = b (1) 

with a set of given complex variables {z^^\ z'''^\ ....z^^^ that have 
physical meaning of energy points. See Sec. 2 for mathematical 
foundation. Since the matrix (z^^^/ — H) is complex symmetric, 
one can solve the equations by the COCG method, indepen- 
dently among the energy points. [8] In these calculations, the 
procedure of matrix-vector multiplications governs the compu- 
tational time. 

For the problem of Eq.(l), a novel Krylov-subspace algo- 
rithm, the shifted COCG method, was constructed [3] [4], in 
which we should solve the equation actually only at one en- 
ergy point (reference system). The solutions of the other energy 
points (shifted systems) can be given without any matrix-vector 
multiplication, which leads to a drastic reduction of computa- 
tional time. The key feature of the shifted COCG method stems 
from the fact that the residual vectors r^^^ = {z^^^ I — H)x{z'^^'>) ~ 
b are collinear among energy points, owing to the theorem of 
collinear residuals. [9] 

Moreover, the shifted COCG method gives another drastic 
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reduction of computational time, when one does not need all 
the elements of the solution vector 07(2:^^^), as in many cases of 
our studies. [3] For example, the inner product 

(2) 

is particularly interested in our cases. The shifted COCG method 
gives an iterative algorithm for the scaler (6, 37(2;^^))), without 
calculating the vector £c(2;*^^^), among the shifted systems. The 
quantity of Eq.(2) is known as 'local density of states' (See 
Sec. 2). It is an energy-resolved electron distribution at a point 
in real space and can be measured experimentally as a bias- 
dependent image of scanning tunneling microscope (See text- 
books of condensed matter physics). 

2 Note on mathematical formulation 

Here a brief note is devoted to mathematical relationship be- 
tween eigen-value equation and shifted linear equation in the 
electronic structure theory. A notation, known as bra-ket no- 
tation in quantum mechanics, is used. See Appendix for the 
details of the notation. 

2.1 Original problem 

Our problem in electronic structure calculation is, convention- 
ally, reduced to an eigen-value problem, an effective Schrodinger 
equation, 

H\Va) = Sa\Va), (3) 

where is a given N x N real-symmetric matrix, called Hamil- 
tonian matrix. Eigen vectors form a complete orthogonal basis 
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set; 



{Va\Va) = 5a(i (4) 

El^a)(^a|=^ (5) 

a 

where / is the unit matrix. 

In actual calculation, the matrix H is sparse. Each basis 
of the matrix and the vectors corresponds to the wavefunction 
localized in real space. Hereafter physical discussion is given 
in the case that the physical system has atoms and only 
one basis is considered for one atom. In short, the i-th basis 
{i = 1,2,3...A^) corresponds to the basis locahzed on the z-th 
atom. Moreover we suppose, for simplicity, that the eigen values 
are not degenerated (ei < £2 < £3 )• 

On the other hand, the linear equation of Eq. ([T]) is rewritten 
in the present notation; 

{z - H)\x,{z)) = \j) (6) 

with a complex valuable z = £ + ^77. The valuable z corresponds 
to the energy with a tiny imaginary part r] [t] ^ +0). 

The purpose within the present calculation procedure is to 
obtain selected elements of the following matrix D; 

D{£) ^6{£-H)=J: \Va)6{£ - £a){Va\ (7) 

a 

or 

Aj(e) = T.{i\Va)5{£ - £a){Va\j)- (8) 

a 

This matrix is called density-of-states (DOS) matrix, since its 
trace 

TT[D{£)]=Y.5{£-£a): (9) 

a 
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is called density of states. The physical meaning of Eq. (Ej) is 
the spectrum of eigen-value distribution. 

In actual numerical calculation, the delta function in Eqs. (tTj) 
and dHj) is replaced by an analytic function, a 'smoothed' delta 
function, since numerical calculation should be free from the sin- 
gularity of the exact delta function ((5(0) = oo). The 'smoothed' 
delta function is defined as 



6^{£) = 



1 



Tm 



TT 



1 



£ + if] 



1 



rj 



TT £^ + rj^ 



(10) 



with a finite positive value of r](> 0). The smoothed delta func- 
tion gives the exact delta function in the limit of 



lim = (5(e). 



The physical meaning of rj is the width of the 'smoothed' delta 
function (5^(e). Hereafter, the DOS matrix is defined as 



D{£) = 5^{£ - H) =J2 \Va)STj{£ - £a){Vc 



(12) 



or 



Aj(e) = E(^l^a)^r,(£ - £a){Va\3)- 



(13) 



It is noteworthy that the smoothed delta function has a finite 
maximum 



^,(0) = - 

T] 

and its integration gives the unity 



/•oo 1 

/ 6J£)d£ = — 
J—oo 7f 



tan 



-1 



£ 



(14) 



£ = 00 



= 1. 



(15) 



e=—oo 
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2.2 Green's function and calculation procedure 

The Green's function is defined as an inverse matrix of 

GW = ^ = E^^^- (16) 

The solution of Eq. (El) gives the Green's function; 

G,,{z)^(^\G{z)\3) = l^\x,[z)) (17) 
and the DOS matrix is given from the Green's function; 

B{e) = --Im [G{£ + ir])] , (18) 

TT 

under the relations of Eq. (1101). In conclusion, the calculation 
procedure, from the Hamiltonian matrix to the DOS matrix, is 
illustrated, as follows; 



H 



PgSd (19) 



2.3 Physical quantities for energy decomposition 

Now we present two quantities, local density of states (LDOS) 
and crystal orbital Hamiltonian population (COHP)[10], as ex- 
amples of important physical quantities that is calculated from 
selected element of the DOS matrix. These quantities appear in 
the decomposition methods of the electronic structure energy. 
The electronic structure energy is defined as 

E = j:^af{Sa). (20) 

a 

where f{£a) is the number of electrons that occupy the wave- 
function with the energy of Sa- In the present case, the case 
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of zero temperature, the function is reduced to a step-function 
form of 

f{e)^e{f,-e) (21) 

with a given value of n (chemical potential). Eqs.dSUj), (El) lead 
us to the expression of the energy with DOS; 

E = r f{E)eY.5{e-Ea)de 

J — OO n 



nOO 



= r f{e)eTT[D{E)]de (22) 

J—oo 

2.3.1 Decomposition with local density of states 
LDOS is defined as diagonal elements 

n,(£) = {i\D{e)\i) = J:\{^K)f5{£ - £a) (23) 

a 

of the DOS matrix and the energy is decomposed into the con- 
tributions of LDOS, {ni{£)}i] 

E = Y.r fi^)^ri,{E)dE. (24) 

• J —OO 

I 

Physical meaning of LDOS is a weighted eigen-value distribu- 
tion; For example, if an eigen vector \va) has a large weight on 
the i-th basis, the local DOS ni(e) has a large peak at the energy 
level of £ = £a- The LDOS ni{£) corresponds to the experimen- 
tal image of the scanning tunneling microscope (See the end of 
Sec. 1). 

2.3.2 Decomposition with crystal orbital Hamiltonian population 

Another decomposition of the energy can be derived with an 
expression of 

/OO 
/(e)Tr[L»(e)//]de. (25) 
-OO 



Eq.([2SD is proved from the first line of Eq.([22I) and the relation 
of 



Tv[D{£)H] 



Tt[6,{£ - H)H] 



a 



a 



a 



(26) 



a 



The last equality is satisfied by the exact delta function (77 
0+). Eq.([25]) gives another decomposition of the energy 



and is called crystalline orbital Hamiltonian population (COHP) 
[10]. The physical meaning of COHP is an energy spectrum of 
electronic wavefunctions, or 'chemical bond', that lie between 
z-th and j-th bases. See the papers [3, 10] for details. 

2.4 Numerical aspects with Krylov subspace theory 

Several numerical aspects are discussed for the calculated quan- 
tities within Krylov subspace theory, such as the shifted COCG 
algorithm. LDOS is focused on, as an example. When we 
solve the shifted linear equation of Eq. (Ej) within the v-ih order 
Krylov subspace 



K^^\H- = span H\j), H'\j), , H^-'\j)} , (29) 





where the matrix C is defined as 




(28) 
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the resultant LDOS should include deviation from the properties 
in the previous sections, properties of the exact solution. 

First, the number of peaks in the LDOS function (^^(e)) is 
equal to the dimension of Krylov subspace, whereas the num- 
ber of peaks in the exact solution, given in Eq.([23j) is N, the 
dimension of the original matrix. A typical behavior is seen 
Fig. 3(a) of Ref.[3], a case with z/ = 30 and = 1024. Here we 
note that the calculation with such a small subspace {v = 30) 
gives satisfactory results in several physical quantities [3] , mainly 
because many physical quantities are defined by a contour inte- 
gral with respect to the energy, as in Eq. (1221). and the informa- 
tion of individual peaks is not essential. 

Second, the calculated function ni(e) in the Krylov subspace 
can be negative (ni(e) < 0), whereas the exact one, given in 
Eq.(l23l), is always positive. In the exact solution of Eq.(l23l), 
peaks in ni(e) are given by the poles of the Green's function 
G{z) {z = £1,^2, ••••) and the function ni{£) contains smoothed 
delta functions of 



Since Eq, is an eigen value of the real-symmetric matrix H and is 
real (Imfe^] = 0), the exact Green's function G{z) has poles only 
on real axis. In the calculation within Krylov subspace, however. 



In that case, the sign of the smoothed delta function can be 
negative; 




(30) 



the poles can be deviated from real axis (eq, = 



w+z4),4)7^o). 




1 



— Im 

TT 



(£-£^^)+^(t/-££^) 
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The above expression indicates that a negative peak appears 
{Sr^{£ — £a) < 0), when the imaginary part of a pole is larger 
than the value oi r] [r] < e'^^). Therefore, we can avoid the 
negative peaks, when we use a larger value of 77, which was 
confirmed among actual numerical calculations with the shifted 
COCG algorithm. 

Finally, we comment on the present discussion from the inter- 
disciplinary viewpoints between physics and mathematics; From 
the mathematical view point, we can say that the above two nu- 
merical aspects disappear, when the dimension of the Krylov 
subspace {v) increases to be enough large. We would like to em- 
phasis that, even if the dimension is rather small, we can obtain 
fruitful quantitative discussion for several physical quantities, as 
discussed above. In other words, mathematics gives a rigorous 
way (iterative solver) to the exact solution and physics gives a 
practical measurement of the convergence criteria. 

A Notation used in Section 2 

In Sec. 2, we use the vector notations of 

l/)^(/i,/2 fuf (32) 

{9\^{9i^92 9m)- (33) 

Particularly, the unit vector of which non-zero component is 
only the i-ih one is denoted as |z); 

|i)^(0,0,....,li,0,0)T. (34) 
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Inner products are described as 



{9\f)^T.9^f^ (35) 

i 

are 

{9\A\f)^T.9^AJfJ (36) 

with a N X N matrix A. The notation of \f){9\ indicates a 
matrix, of which a component is given as 

{\f){9\).j = f^9j. (37) 

These notations are used in quantum mechanics and called 'bra- 
kef notation. We should say, however, that the above notations 
are slightly different of the original 'bra-kef notations. For ex- 
ample, the original notation of {g\ is {g\ <^ {9i,92 9m)- The 

reason of the difference comes from the fact that the standard 
quantum mechanics is given within linear algebra with Hermi- 
tian matrices but the present formulation is not. 
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